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ABSTRACT 

We  consider  using  data  acquired  from  a  micro-lens  array  through  which  multiple  images  of  the  full 
field-of-view  of  an  astronomical  target  are  formed  to  attempt  to  reconstruct  the  3-D  wave  front  for  the 
observations.  This  opens  the  door  for  both  a  beacon-less  wave  front  sensor  and  imaging  of  fields-of- 
view  substantially  larger  than  the  isoplanatic  angle.  The  reconstruction  problem  can  be  modeled  as  a 
large-scale  linear  inverse  problem,  but  standard  algorithms  used  for  3-D  computed  tomography  (CT) 
reconstruction  cannot  be  applied  because  measured  data  are  only  taken  from  limited  angular  range, 
leaving  entire  regions  of  the  frequency  space  un-sampled.  However,  we  show  that  there  is  substantial 
structure  in  the  mathematical  model  that  can  be  exploited  to  obtain  a  robust  algorithm  that  is  amenable 
to  efficient  implementations. 


1.  INTRODUCTION 

In  many  imaging  situations,  such  as  when  ground  based  telescopes  are  used  to  observe  objects  in  space, 
the  observed  image  is  degraded  by  blurring  and  noise.  Although  the  blurring  can  be  partially  removed 
through  sophisticated  (and  expensive)  imaging  devices,  such  as  adaptive  optics  telescopes,  computational 
postprocessing  techniques  are  also  often  needed  to  further  improve  the  resolution  of  the  image.  This 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

SEP  2014 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

Fast  Tomographic  Reconstruction  of  Atmospheric  Turbulence  from 
Micro- Lens  Imagery 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Emory  University, Mathematics  and  Computer  Science, Atlanta, GA, 30322 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2014  to  00-00-2014 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

In  the  Advanced  Maui  Optical  and  Space  Surveillance  Technologies  (AMOS)  Conference,  9-12  Sep  2014, 
Maui,  HI. 

14.  ABSTRACT 

We  consider  using  data  acquired  from  a  micro-lens  array  through  which  multiple  images  of  the  full 
eld-of-view  of  an  astronomical  target  are  formed  to  attempt  to  reconstruct  the  3-D  wave  front  for  the 
observations.  This  opens  the  door  for  both  a  beacon-less  wave  front  sensor  and  imaging  of  elds-of-  view 
substantially  larger  than  the  isoplanatic  angle.  The  reconstruction  problem  can  be  modeled  as  a  large-scale 
linear  inverse  problem,  but  standard  algorithms  used  for  3-D  computed  tomography  (CT)  reconstruction 
cannot  be  applied  because  measured  data  are  only  taken  from  limited  angular  range  leaving  entire  regions 
of  the  frequency  space  un-sampled.  However,  we  show  that  there  is  substantial  structure  in  the 
mathematical  model  that  can  be  exploited  to  obtain  a  robust  algorithm  that  is  amenable  to  e  cient 
implementations. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

9 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


computational  postprocessing,  which  is  referred  to  as  deblurring,  restoration,  or  deconvolution  [1,  5,  12, 
10,  13,  21],  requires  solving  an  ill-posed  inverse  problem 


g(x,y) 


ty-  V  )/(£>  y)d^drj  +  e(x,  y ) , 


(1) 


where  /  is  the  true  object,  g  is  the  observed  image,  and  e  is  additive  noise.  The  kernel  function  k  models 
the  blurring  operation,  and  is  called  the  point  spread  function  (PSF).  Using  a  standard  Fourier  optics 
model  for  atmospheric  turbulence  [18],  the  PSF  can  be  expressed  in  terms  of  the  wavefront  phase  error, 
</>,  of  the  light  that  reaches  the  telescope  mirror, 


k{x,y) 


T~x  ^P{x,y)e^x’y)^ 


(2) 


where  l  =  yf—1,  V  is  a  characteristic  function  that  models  the  shape  of  the  telescope  aperture  (e.g.,  a 
circle  or  annulus),  and  T~l  is  the  2-dimensional  inverse  Fourier  transform.  In  perfect  seeing  conditions, 
cf>  =  0,  and  in  good  seeing  conditions  (low  turbulence),  </>  is  a  smooth  function. 

The  effectiveness  of  image  restoration  algorithms  depends  strongly  on  how  well  one  can  approximate 
the  PSF,  or,  equivalently,  how  well  one  can  estimate  the  phase,  f>.  In  some  cases,  a  wavefront  sensor 
(WFS)  can  be  used  to  measure  gradients  of  the  wavefront,  from  which  the  phase  can  be  reconstructed. 
A  WFS  is  standard  technology  in  adaptive  optics  systems,  and  many  papers  have  been  written  about 
efficiently  reconstructing  the  wavefront  from  the  gradient  measurements;  see,  for  example  [3,  11].  The 
phase  reconstruction  problem,  from  measured  gradients,  depends  on  the  sensor  geometry  [6,  11]  but  it 
essentially  amounts  to  solving  a  discrete  finite  difference  problem  [8,  17,  4].  Probably  the  most  commonly 
used  WFS  is  the  Shack-Hartmann,  found  in  optical  systems  across  a  very  broad  range  of  applications  in 
astronomy,  medical  imaging,  space  situational  awareness,  directed  energy  weapons,  and  secure  commu¬ 
nications  [22,  9,  7].  But  the  Shack-Hartmann  requires  a  well-defined  and  preferably  unresolved  source  on 
the  far  side  of  the  turbulence  to  act  as  a  beacon.  In  space  imaging  applications,  this  is  either  a  natural 
guide  star  or  a  laser  guide  star. 

The  wavefront  phase  gradient  information  measured  by  the  WFS  is  an  integration  of  the  turbulence 
effects  along  a  line  of  sight  through  the  3-dinrensional  atmosphere.  If  imaging  over  a  narrow  field  of 
view  (FOV),  then  it  is  generally  not  necessary  to  model  the  3D  turbulence  profile;  a  2D  integrated  phase 
provides  sufficient  information  to  construct  a  spatially  invariant  PSF.  However,  when  imaging  over  a  wide 
FOV,  modeling  the  3D  atmospheric  structure  is  important  because  a  spatially  invariant  assumption  on 
the  PSF  may  not  be  appropriate,  especially  in  extreme  turbulence  situations. 

Atmospheric  tomography  techniques  have  been  developed  to  reconstruct  a  pseudo-3D  representation 
of  atmospheric  turbulence.  We  use  the  phrase  pseudo-SD  because  it  is  not  feasible  to  obtain  a  full 
3D  reconstruction  of  the  atmosphere.  Instead,  the  aim  of  atmospheric  tomography  is  to  reconstruct 
information  at  a  discrete  set  of  dominant  layers  (slices)  of  atmospheric  turbulence.  In  order  to  do  this, 
it  is  necessary  to  use  multiple  measurements,  where,  loosely  speaking,  each  measurement  is  obtained 
by  viewing  the  turbulence  layers  from  a  different  perspective.  The  typical,  and  perhaps  most  obvious, 
approach  to  do  this,  is  by  collecting  WFS  measurements  of  guide  stars  distributed  throughout  the  FOV, 
and  use  standard  limited  angle  tomography  techniques  to  reconstruction  the  turbulent  layers  [2,  15,  16, 
19,  20,  23,  14].  This  is  the  approach  taken  in  the  Gemini  multi-conjugate  adaptive  optics  system  [14], 
which  employs  five  sodium  laser  guide  stars  to  correct  a  field  of  view  of  about  an  arcminute. 

Because  it  may  be  impractical  or  undesirable  to  obtain  multiple  observations  of  guide  stars,  we  take 
an  alternative  approach.  Specifically,  we  consider  using  WFS  data  acquired  from  a  micro-lens  array  in 
which  the  full  FOV  of  an  astronomical  target  is  imaged  through  each  micro-lens.  Gradients  of  the  local 
wavefront  phase  aberration  over  each  subaperture  are  used  to  reconstruct  phase  gradients  at  a  discrete 
set  of  turbulence  layers.  With  known  phase  gradients,  it  is  then  straightforward  to  reconstruct  phases  at 
each  layer. 


We  show  that  the  problem  of  reconstructing  phase  gradients  at  each  layer,  given  phase  gradients 
in  each  subaperture,  can  be  modeled  as  a  large-scale  linear  inverse  problem.  The  problem  is  severely 
ill-conditioned,  and  can  have  distinct  null  space  components.  However,  in  certain  situations,  and  with 
appropriate  regularization,  we  show  that  phases  can  be  reconstructed  at  multiple  layers. 

The  rest  of  this  paper  is  organized  as  follows.  In  Section  2  we  describe  the  problem  setup,  discuss 
the  properties  of  the  mathematical  model,  and  our  approach  to  solve  the  resulting  inverse  problem.  In 
Section  3  we  show  the  results  of  some  experiments  on  simulated  data  to  illustrate  the  effectiveness  of  our 
approach.  Some  concluding  remarks  are  given  in  Section  4. 

2.  PROBLEM  SETUP 

We  consider  a  situation  where  the  WFS  divides  the  telescope  pupil  into  a  set  of  discrete  lenses;  see 
Figure  1  for  an  illustration  of  5  lenses  (one  should  really  think  of  this  as  a  discrete  grid  of  5  x  5  lenses, 
but  the  one-dimensional  grid,  or  side  view  of  a  two-dimensional  grid,  shown  in  Figure  1,  simplifies  the 
notation).  Each  subaperture  captures  light  reflected  off  of  objects  in  the  FOV,  but  from  slightly  different 
perspectives.  In  particular,  as  light  travels  through  the  atmosphere,  it  is  affected  by  the  dominant  layers 
of  turbulence  in  a  cone  region  defined  by  the  diameter  of  the  FOV  and  the  diameter  of  the  subaperture. 
These  cones  are  illustrated  in  Figure  1  with  different  colored  lines;  for  example,  red  lines  depict  the  cone 
region  corresponding  to  the  first  (far  left)  subaperture. 

Equation  (2)  can  be  used  to  model  local  blurring  in  the  images  captured  by  each  subaperture,  but 
because  of  the  slightly  different  perspectives  of  the  subapertures,  each  local  blurring  is  defined  by  a 
different  phase.  Specifically,  the  phase  corresponding  to  the  i- th  subaperture  is  defined  by  integrating 
through  the  atmospheric  layers  within  the  cone  region  from  the  *-th  subaperture  to  the  FOV.  Note 
that  the  WFS  actually  measures  gradients,  and  not  phases  directly,  but  since  differentiation  is  a  linear 
operation,  we  can  assume  the  WFS  measured  gradients  are  an  integration  of  gradients  in  the  same  cone 
regions. 

To  describe  a  mathematical  model  that  relates  the  measured  phase  gradients  to  phase  gradients  at  the 
various  layers  of  turbulence,  let  (j>x{i )  denote  measured  phase  gradients  in  the  x,  or  horizontal  direction, 
and  4>v{i)  denote  measured  phase  gradients  in  the  y,  or  vertical  direction,  in  the  *-th  subaperture.  We 
will  assume  an  na  x  ns  square  grid  of  subapertures,  and  define  Ns  =  n2sl  and  thus  i  =  1,  2, . . . ,  Ns. 

Suppose  further  that  we  partition  the  unknown  phase  gradients  at  the  various  layers  of  turbulence 
into  an  np  x  np  grid,  and  denote  these  partitioned  regions  of  phase  gradients  as  <j>x(j ,  £)  and  <j)y(j,  £),  where 
j  =  1,  2, . . . ,  Np  =  np,  and  l  =  1,2,...  Nl,  with  Nl  defining  the  number  of  layers  of  turbulence  that  we 
want  to  reconstruct.  Figure  2  illustrates  our  notation,  which  again  shows  a  side  view  of  a  two-dimensional 
grid  with  Np  =  4,  Nl  =  3,  and  Ns  =  25. 

Using  this  notation,  a  mathematical  model  relating  the  measured  phase  gradients  to  the  unknown 
phase  gradients  at  each  atmospheric  layer  can  be  given  by 

nl  np 

=  'Y^2,w{i,j,£)<f>x{j,£)  +£x(i) 

t-i  j- i 

V  (3) 

Nl  Np 

e=i  j= i 

where  w(i,j,£)  is  a  weight  associated  with  the  contribution  of  to  and  £x(i)  and  £y(i)  represent 

unknown  noise  and  other  measurement  errors.  The  weight  w(i,j,£)  is  the  fraction  of  the  region  in  which 
<f)(j,  £)  is  defined,  which  is  contained  within  the  i- th  cone.  For  example,  the  situation  depicted  in  Figure  2 
shows  that  w( 2,  2, 1)  =  0.  In  general,  it  is  expected  that  for  a  given  subaperture,  i,  most  of  the  weights 
w(i,j,£)  will  be  zero. 
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Figure  1:  This  figure  illustrates  a  situation  in  which  there  are  three  dominant  atmospheric  layers  of 
turbulence  (one  at  the  ground,  and  two  at  higher  altitudes),  and  5  lenslets  in  the  aperture.  In  reality, 
there  should  be  a  grid  of  lenslets  (e.g.,  5x5),  but  to  simplify  the  illustration,  we  show  only  a  one- 
dimensional  side  view  of  the  aperture. 


Figure  2:  This  figure  illustrates  discretization  of  the  phase  gradients. 


Note  that  the  equations  in  (3)  can  be  written  in  matrix-vector  form, 

4>x  =  Wcj)x  +  ex  and  <j>y  =  W<j)y  +  ey,  (4) 

where  W  is  an  Ns  x  NpN^  matrix  whose  coefficients  are  given  by  w(i,j,£).  It  is  important  to  observe 
that  to  avoid  an  explicitly  under-determined  system,  we  should  require  that  the  partitioning  of  the  phase 
gradients  on  each  layer  satisfy  Np  <  Ns/Nl ■  Even  with  such  a  restriction,  the  matrix  W  can  be  highly 
ill-conditioned,  with  a  possibly  nontrivial  null  space;  this  is  illustrated  in  Section  3.  Thus,  solutions  can 
be  very  sensitive  to  noise  and  other  data  errors,  and  it  is  necessary  to  incorporate  regularization  in  the 
reconstruction  methods.  In  this  paper,  we  use  standard  Tikhonov  regularization;  that  is,  we  solve  the 
least  squares  problems: 


min 

1 fix 


w 

axI 


X 

o 


and  min 

0 y 


w 

OiyI 


4>y 


(5) 


We  remark  that  although  the  two  problems  in  equation  (5)  have  the  same  coefficient  matrix,  the  fact 
that  they  have  different  data  vectors  means  that  the  “optimal”  level  of  regularization  (defined  by  scalars 
ax  and  ay)  might  be  different.  However,  currently  we  are  using  the  same  level  of  regularization  for  both 
problems,  and  we  leave  for  future  work  a  more  complete  analysis  of  the  two  problems,  as  well  as  the 
design  of  optimal  regularization  methods. 


3.  NUMERICAL  EXPERIMENTS 

In  this  section  we  provide  experimental  results  on  simulated  data  to  test  the  potential  of  our  proposed 
approach  to  atmospheric  tomography.  In  each  simulation,  we  assume  that  the  FOV  is  at  330  km  (the 
altitude  of,  for  example,  the  International  Space  Station),  coverage  area  of  40  arcsecs,  two  dominant 
layers  of  turbulence  (Nl  =  2),  and  an  aperture  diameter  of  3.6  nr. 

We  use  the  forward  models  given  in  equation  (4),  where  we  first  use  ex  =  ey  =  0,  and  then  take  these 
terms  to  be  1%  Gaussian  white  noise;  that  is,  ex  and  ey  are  vectors  with  normally  distributed  random 
entries,  mean  0,  standard  deviation  1,  and  scaled  to  that 

ll^'lh  _  ll^lb  _  „  m 
\\W<t>X\\2  \\W<fiy\\2  ■  ■ 

Note  that  even  without  additional  noise,  the  data  contains  downsampling  approximation  errors,  which 
can  have  a  deleterious  effect  on  the  computed  reconstructions. 


We  consider  two  possible  configurations: 


ns  =  32,  np  =  16  and  ns  =  64,  np  =  32  . 

The  first  configuration,  where  we  assume  the  aperture  is  partitioned  into  32  x  32  subapertures,  reflects 
current  capabilities  of  a  typical  WFS. 

Figure  3  shows  the  singular  values  of  the  reconstruction  matrix,  W,  for  both  configurations.  From 
these  plots  we  see  that  the  matrix  has  a  one-dimensional  (numerical)  null  space,  and  except  for  the 
one  tiny  singular  value,  the  remaining  singular  values  are  relatively  large.  Therefore,  we  expect  that  a 
reasonable  reconstruction  of  the  layers  should  be  possible,  provided  we  incorporate  a  modest  amount  of 
regularization,  and  provided  that  the  null  space  component  of  the  solution  does  not  contribute  significant 
information  to  the  reconstructed  phases. 

Reconstructions  for  the  configuration  with  ns  =  32  and  np  =  16,  using  both  noise  free  and  noisy  data, 
are  shown  in  Figure  4  .  Reconstructions  for  the  configuration  with  ns  =  64  and  np  =  32,  using  both  noise 
free  and  noisy  data,  are  shown  in  Figure  5  .  For  both  configurations  we  used  regularization  parameters 

=  Oiy  =  10-°  for  the  noise  free  data,  and  ax  =  ay  =  10-2  for  the  noisy  data. 


Figure  3:  This  figure  shows  the  singular  values  for  the  matrix  W;  the  left  plot  shows  singular  values 
for  the  case  ns  =  32  and  np  =  16,  and  the  right  figures  shows  singular  values  for  the  case  ns  =  64  and 
np  =  32. 

These  results  illustrate  that  our  model  and  computational  approach  to  solve  the  atmospheric  to¬ 
mography  problem,  while  not  perfect  at  reconstructing  atmospheric  layers,  still  provides  relatively  good 
reconstructions  considering  we  use  only  one  set  of  collected  data  of  a  scene  of  interest.  In  particular,  we 
do  not  do  not  collect  multiple  images,  and  we  do  not  to  use  laser  guide  stars.  It  appears  that  the  null 
space  component  is  important,  and  further  work  is  needed  to  understand  its  significance,  and  how  to 
incorporate  its  information  into  the  solution. 

4.  CONCLUDING  REMARKS 

Previously  proposed  techniques  for  atmospheric  tomography  typically  use  a  distribution  of  guide  stars  to 
obtain  information  from  which  to  obtain  a  pseudo-3D  reconstruction  of  atmospheric  turbulence  layers. 
Because  this  may  not  always  be  practical,  or  feasible,  we  propose  an  alternative  approach  using  a  WFS 
that  is  partitioned  into  a  grid  of  subapertures.  Our  approach  requires  the  collection  of  only  one  set  of 
data  of  a  scene  of  interest,  and  does  not  require  the  use  of  any  laser  guide  stars.  We  provided  a  math¬ 
ematical  and  computational  framework  to  use  the  phase  gradients  of  the  local  blur  in  each  subaperture 


Figure  4:  Comparison  of  true  phase  with  reconstructed  phases  for  the  case  ns  =  32  and  np  =  16.  The 
top  row  is  the  phase  the  layer  at  0  km  and  the  bottom  row  is  the  layer  at  5  km.  The  left  column  is  the 
true  phase,  the  middle  column  is  the  reconstructed  phase  using  noise  free  data,  and  the  right  column  is 
the  reconstruction  using  noisy  data. 


Figure  5:  Comparison  of  true  phase  with  reconstructed  phases  for  the  case  ns  =  64  and  np  =  32.  The 
top  row  is  the  phase  the  layer  at  0  km  and  the  bottom  row  is  the  layer  at  5  km.  The  left  column  is  the 
true  phase,  the  middle  column  is  the  reconstructed  phase  using  noise  free  data,  and  the  right  column  is 
the  reconstruction  using  noisy  data. 


to  reconstruct  phase  gradients  at  a  discrete  set  of  turbulence  layers,  which  is  then  followed  by  a  phase 
reconstruction  on  each  layer.  The  numerical  experiments  presented  in  this  paper  illustrate  that  our  model 
can  be  an  effective  approach  to  solve  the  atmospheric  tomography  problem.  There  are  still  many  open 
problems  that  we  intend  to  address,  including  an  analysis  of  the  reconstruction  matrix,  W ,  and  more 
specifically,  how  to  better  regularize  the  null  space  components. 
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